Some topics in the kinetics of protein aggregation 
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Abstract 



Preliminary results are presented for the kinetics of phase sepa- 
ration in three distinct models of protein aggregation. The first is 
a model of the formation of spherical microcrystals of insulin via an 
initial formation of fractal clusters of insulin. The results of our Brow- 
nian dynamics study of this model are in qualitative agreement with a 
r^ recent experimental study [1 of microcrystal formation from aqueous 

0^ mixtures of insulin. A second work involves a theory for the formation 

of metastable bundles of sickle hemoglobin from fibers, based on a re- 

^^ cent generic theory of bundle formation [2 . We also discuss a model 

^ for the microscopic formation of these fibers. Finally we discuss pre- 

^N liminary results for the kinetics of cluster formation for a six patch 

model of protein crystallization. 



Keywords: protein aggregation, brownian dynamics, sickle hemoglobin, 
patch model 

1 Introduction 



S^ In the spirit of this Workshop, we present some preliminary results for several 



problems that involve the kinetics of protein aggregation. This subject is an 
active field of research that includes studies of protein crystallization as well 
as certain biomedical problems such as human cataract formation, sickle 
cell anemia, Alzheimer's disease and various problems in drug delivery. The 
specific problems discussed below have in common that they involve the 
kinetics of phase separating systems. A review of recent developments in 
this field is contained in 131. 



2 A Model of Microcrystal Formation in Insulin 
Solutions 

A standard method of preparing microcrystals of insulin for use in drug de- 
liver is to precipitate insulin from aqueous solutions using zinc salts. This 
is the same technique as is used in many kinds of protein crystallization. 
However, in 2003 Bromberg et al [I] used an alternative method in which 
the aqueous solution was prepared near the isoelectric point of insulin, in 
order to minimize Coulomb interactions. Rather than using a salt to pre- 
cipitate insulin from solution, they used polyethylene glycol (PEG), which 
is known to induce attractive interactions between biomolecules via a deple- 
tion attraction. The authors first quenched the system to a low temperature 
and found that the insulin particles precipitated from solution in the form 
of a fractal network. The fractal nature was established by continuous- angle 
small angle light scattering, which showed a fractal dimension of 1.8. This 
dimension is characteristic of the case of diffusion-limited cluster-cluster ag- 
gregation. They then stirred and subsequently diluted the sample to break 
up the fractal network into a relatively monodisperse set of microcrystals. 
This promises to be an alternative method of microcrystal preparation in 
the drug delivery of insulin. 

Our group has recently carried out a Brownian dynamics simulation of 
a model of this experiment, to see if we are able to capture the qualitative 
features of the experiment. Our Hamiltonian consists of a hard core inter- 
action between the insulin molecules, together with a short-range attractive 
interaction given by the Asakura-Oosawa model of depletion attraction. The 
latter is a semi-quantitative approximation of the depletion interaction in- 
duced by the PEG polymers. In our Brownian dynamics (BD) simulations 
, we consider a three-dimensional system of size L — 64cr in units of the 
insulin molecular diameter a. All other length scales are measured in units 
of a as well. We consider the case of low volume fraction / = 0.02 for a 
system of 10,013 insulin molecules. The equations of motion for the BD 
simulation are 

F, = -vt/,-rf, + w>,(t) (1) 

where F is the friction coefficient and W^, the random force acting on each in- 
sulin particle i, is a Gaussian white noise satisfying a fluctuation-dissipation 
relation. Hydrodynamic interactions, including lubrication forces are ig- 
nored in the simulation. The potential U acting upon each insulin monomer 
has a twofold contribution: the two-body depletion potential of the Asakura- 



Oosawa-Vrij (Uao) plus a repulsive hard-core-like interaction (Uhc) given 
by the following expressions: 



u{nj) = UAoiuj) + Uncinj) (2) 



where 



and is zero for r^j > (1 + ^). The hard core potential is given by 






= ^r • (4) 



u 



In Eq. 3, ^ is the size-ratio between a polymer coil and a colloidal 
particle which controls the range of the depletion interaction in the Asakura- 
Oosawa-Vrij model and (j)p is the polymer volume fraction which controls the 
strength of the interaction. Our simulations are for ^ = 0.1. In the hardcore- 
like repulsive interaction given by Eq. 4, we have set a = 36. Values of 
a < 36 have been reported to lead to anomalies when a mimic of the hard- 
core potential is required in the potential [18,33]. The total pair potential 
U = Uao + Uhc has a minimum value (Um) that depends on ^ and 0^. 
In what follows, we will often characterize the strength of the potential in 
terms of the absolute value of the minimum potential depth, Um = |^mm|- 
We choose F = 0.5 and a time step At = 0.005 in reduced time units of 
a{m/kBT)^/'^, with m = 1. For this choice of F, particle motion is diffusive 
for t » p, i.e.t ^ 2 in our units. Periodic boundary conditions are used to 
minimize wall effects. All simulations start from a random initial monomer 
conformation and the results for the kinetics are averaged over several runs. 

2.1 Results 

We summarize here some of the results of our simulation. We first quench 
the system deep into the two phase gas-solid region (with \Um\ = lO.O/c^T) 
and study the kinetics of the resultant cluster formation. Figure [l] shows 
the morphology of the system for various time steps. It appears that the 
system initially phase separates through the formation of fractal clusters. 
We have verified this in more detail in various ways; Figure [2] shows the 
behavior of the number of clusters, Nc and the radius of gyration, Rg as a 
function of time, in a log-log plot. These have slopes of —1 and 0.55, respec- 
tively, in the early stages of development, in agreement with the diffusion 




Figure 1: Morphology of cluster formation (^ = 0.1, / = 0.02) for a deep 
quench {\Um\ = lO.O/c^T) into the two phase gas-solid region at various 
times: (a) t = 0, (b) t = 50, (c) t = 1, 000, and (d) t = 10, 000. 



limited cluster-cluster aggregation theory (DLCA). We have also studied 
the structure factor of the system at various times (see Figure [s]) and found 
that this behaves like S{q) ^ q~^^ ^ where Df = 1.8, for small q, consis- 
tent with fractals with a dimension Df = 1.8. The large q behavior is 
consistent with Porod's law, which characterizes the scattering from com- 
pact clusters, namely S{q) ^ q~^^^^\ d = 3. Thus the clusters are hybrid 
fractals with short-range crystalline order and long-range fractal morphol- 
ogy. To simulate the experimental situation in which the insulin mixture 
is stirred and diluted after the formation of fractals, we have taken out the 
largest cluster, put it into another simulation box and "heated" the system 
to \Um\ = 2.88A:^r. Its subsequent time evolution is shown in Figure [4J 
As can be seen this process results in a break-up of the fractal cluster into 
spherical aggregates, as shown, say, in panel c) of that figure. We then 
"cooled" the system in b) to \Um\ = 3.72kBT to further stabilize this dis- 
tribution of droplets. By this process, we have been able to show that we 
can reproduce the essential features of the experimental study by Bromberg 
et al [Ij via the Asakura-Oosawa model for depletion attractions induced 
by PEG. A more quantitative theory of their experiment would require the 
inclusion of additional forces that would lead to a temperature dependent 
behavior of the system, as seen experimentally. 




Figure 2: Plot of (a) number of clusters, Nc, and (b) radius of gyration, 
Rg, as a function of time in a log-log plot, for a lOfc^T, / = 0.02, ^ = 0.1 
system; this shows slopes of —1 and 0.55, respectively, in the early stages, 
in agreement with DLCA. 




Figure 3: (a) Structure factors and (b) log-log plot of structure factors at 
several times for \Um\ = lO.Ofc^T. Dashed line indicates fractal clusters with 
S{q) ^ q~^^ , where Df = 1.8. The dotted line indicates the Porod regime 
S{q) ^ q~^^^^\ d = 3. Thus the clusters are hybrid fractals with short-range 
crystalline order and long-range fractal morphology. 




Figure 4: (a) Largest cluster from a deep quench \Um\ = lO.Ofc^T, ^ = 0.1, 
t = 10, 000. Cluster morphology at t = 5, 000 for break-up of the fractal 
clusters and formation of spherical aggregates, (b) Heating up the entire 
system (a) to \Um\ = 2.88A:^r for t = 5,000; (c) keeping the entire system 
in (b) at \Um\ = 2.88/c^r for another t = 5,000 and (d) cooling the entire 
system in (b) to \Um\ = 3.72kBT for t = 5, 000. 



3 A Theory for the Finite Bundle Size of Sickle 
Hemoglobin Molecules 

Aggregates, or bundles, of twisted protein fibers, such as sickle hemoglobin 
and actin, are important examples of biopolymers in which elastic interac- 
tions play a crucial role in determining the (metastable) bundle radii. In 
one recent paper [3| Turner et al. proposed a model for stabilizing approx- 
imately 20 nm diameter bundles of sickle hemoglobin (HbS) fibers. They 
constructed the free energy per unit volume, G, needed to create a fiber 
bundle, where G = F — ^, using continuum elasticity theory. Here F is the 
distortion free energy per unit volume of a bundle of radius R and pitch 
length A and '0 is a positive Lagrange multiplier that controls the pitch 
length. From G they predicted the physical properties of the fiber bundle, 
such as the equilibrium (metastable) bundle radius Re, where in the latter 
case they minimized G with respect to R. However, we believe their analysis 
is incorrect for two reasons, the first being the use of the free energy density, 
G, rather than the total free energy E?LG, to determine Re. The second 
is their omission of the binding energy between fibers, which in classical 
nucleation theory of spherical droplets corresponds to the driving force for 
nucleation. We present a corrected version of their analysis below. Our 
approach is the same as that of Grason and Bruinsma L2j, who determined 
the critical bundle size for aggregates of filamentous actin. 
According to classical homogeneous nucleation theory [2l|5l|6j, the critical 
"droplet" size corresponds to the maximum of the total free energy R^LG, 
which is significantly different from the maximum of the free energy density 
G. A simple example is the nucleation of a spherical droplet [6], in which the 
total free energy is given by 47ri?^7 — ^i?^e, where 7 is the surface tension 
and e is the free energy density difference between the metastable and stable 
phases. G{R) — ^ —^- The critical radius Re is determined by maximizing 
this total free energy (not the free energy per unit volume) with respect to 
R. The analogous argument for the heterogeneous nucleation of the fiber 
bundle involves calculating the total free energy involved in creating this 
bundle from an aggregate of fibers of (fixed) length L. 

The grand potential of a twisted fiber ^q as a function of pitch A and 
radius R, includes the contributions from the surface tension, extension or 



f = 10/m"\7 = 5.4///m"' 




yi^o ■ 7^;o 



Figure 5: Plot of the free energy ri(i?; A*) per unit length as a function of 
the fiber radius, R, using the experimental values for HbS given in the text 
for e > 190J7Ti~^ (black), e = 38 Jm~^ (red), and e = IOJtti"^ (green). One 
local minimum occurs at i? = llnm which corresponds to the (metastable) 
equilibrium radius of HbS when e = 38J7ti~^,7 = 5.8/xJ7n~^. 



compression, bending, twisting, binding and chemical potential: 






?2 16 ^ 96 ^ ^ d2^ 



f^ = 7vL{2-fR + ER' ^^ ^^ -^-R'e) (5) 

where L is the fiber length, E the extensional modulus, a the radius of a 
protofilament, and fi the chemical potential of protofilaments. ^ is related 
to the twisting stiffness %A\' Equation [s] contains an additional term — i?^e 
due to the aggregation energy [5l [6] between fibers that is not present in 
Turner et al. [4J. In the limit of L ^ oc, the grand potential is dominated 
by ri. Therefore the equilibrium pitch is determined by gf^^A) Ia=a* — 0? 
which reduces Q to 



27i?-— 7^— ^— ^-i?2 (6) 



ttL ' 27r4/3£;i/3 (6a2 + i?2)l/3 

)erimental 
^-' H [7 
e > 190 Jm-3 (Fig. |5) 



Using experimental values for HbS of a = 4n77i, E — 51MPa, '^ = 
3.5 X 10~^ Jm~^ [U ^, we find that f^(i?; A*) has just a single peak for 

R — {) and R ^ oo correspond to the phases of 
the dispersed protofilaments and stable crystal structures, respectively. As 
e decreases below this, a local minimum develops in Q(R; A*) whose position 
depends on e and 7. The minimum critical bundle size Re occurs under the 
condition that ri(i?; A*)|i^^ = 0, — ^ — -\r^ = 0. Combining the estimate 
e ^ 38 Jttt."^ for HbS [7J, this yields a value of Re = llnm and 7 = 5.8fiJm~^ 
( Fig. [5]), which are consistent with experimental observations [HE]. A 
further reduction in e leads to a decreasing value of ^{Rc] A*) (Fig. [5]). We 
also note that the torsional rigidity obtained by Turner et al is the same 
in our calculation, because we use their approximation for the elastic free 
energy; this value for the rigidity is in agreement with experimental values. 
Finally, there always is an energy barrier for the transition from dispersed 
protofilaments to the metastable bundle phase, which is incorrectly predicted 
as a spontaneous process in reference [4j. 

4 Preliminary Results for a Brownian Dynamics 
Simulation of Bundle Formation 

To understand and compare to experimental observations and our theo- 
retical predictions of HbS, we are carrying out Brownian dynamics (BD) 



simulations of bundle formation based on microscopic interactions. In our 
simulation, the chiral filaments are described by the helical wormlike chain 
model [8j, in which the bending and twisting energy are incorporated into 
bead-spring polymers. The potential energy U acting upon each monomer 
has three contributions: the elastic energy associated with a single chain, 
the repulsive energy due to the excluded volume and the highly anisotropic 
short-range attractive energy between chains. 

^ — ^ chain ~r t-^rep i ^ c—c 



1 N N N 

Uchain = y X^Ki-1 - lof + y X](^i - Ui-lf + y ^(t^ 
i=l i=l i=l 



ro? 



N 

where kg is the spring constant, N the polymer length, Vi^i-i the distance 
between two adjacent monomers ( ith and (i — l)th), /^^ the bending stiffness. 
Hit the torsional stiffness, Ui the tangent vector and r^ the torsional angle on 
ith monomer, r^j is the distance between two monomers, a the diameter of 
each monomer and e^j the pair well depth. The anisotropic attraction is 
modulated by the patchy model [9J: 

1 ^ 

Uc-c = 2 ^ Uattr{rij)Vang{rij,^iMj) (8) 

where Uattr is a Yukawa potential 

C _^^ eM-Z{n,/a-l)) if >^ 

Uattr = ^ '- .^ '' - (9) 

[0 it Tij < a. 

Here A is the energy strength and Z characterizes the range of attraction. 

q2 q2 

Vanginj^^i.^j) = exp{-^)exp{-^) (10) 



where O^^ij is the angle between patch k on the ith monomer and the in- 
terparticle displacement fij. The particular pair of patches chosen (of the 
eight possible patches) is that which minimizes the magnitude of the angles 
Ok^ij and Oi^ij respectively. The parameter ^o is the standard deviation of the 
Gaussian distribution. As an example, we apply the model for two simple 
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(a) (b) 

Figure 6: (a) a single straight chain with twenty monomers at t = 0, (b) a 
hehcal structure is formed at t = 10,000. 



cases: a single chiral chain and two binding chiral chains. In the first case, 
a single straight chain is initially introduced in the system (Fig. |6(a)[). As 
time evolves, a chain twists into a helical structure, shown in Fig. |6(b)| 

In the second case of two binding chains, four patches are symmetri- 
cally arranged on the equator (perpendicular to the polymer bond) of each 
monomer. Two chains, one straight and one helical, are introduced, in which 
patches are aligned so that the patchy attraction is imposed at t = (Fig. 



7(a)). Afterwards, two chains start intertwining and a configuration at t = 
100 is shown in Fig. |7(b)[ At t = 5000, two chains form a double helical 
structure (Fig. |7(c) ). 



5 Kinetics of Cluster Formation in a 6-Patch Model 
of Protein Crystallization 

Many proteins are globular in shape but they have nonuniformly distributed 
surface charges that yield highly anisotropic interactions. In a recent work, 
the process of self-assembly in protein crystallization has been studied ^lOj . 
In that work, they model spherical proteins using a patch model. The 
interaction between proteins takes into account an isotropic and a highly 
directional interaction. The isotropic part was modeled as a square well 
of depth 6 within a range of interparticle separation between a and Aa, 
(J being the diameter of the particles and A = 1.15. For the anisotropic 
part, six patches uniformly distributed on the surface of each particle was 
considered. Then, the intrapatch interaction potential has a radial and an 
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(a) 



(b) 




(c) 

Figure 7: (a) two bound chains with four equatorial patches (b) a configu- 
ration at t = 100, (c) at t = 5000, a double helical structure is formed. 



angular dependence: 



(11) 



where u^^ is defined in the same way as the isotropic part, but in this case 
with a deeper well depth Cp = 5e and narrow well width Xp = 1.05. The 
angular part was defined as: 



J \^Hi ^^j J 



1 ei,ej<6 

otherwise 



9i being the angle between the patch orientation, determined by a unit 
normal vector centered in the patch, and the line that connect the center of 
mass of the particles. S = 0.259 determines the patch size. 
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In our work, we are interested in the dynamics of the nucleation and 
growth process of a six-patch model of globular protein using Brownian 
dynamics simulations. In particular we study the kinetics of the process by 
which chains form and develop local crystal order that then optimizes the 
crystal nucleation process. 



Figure 8: Sketch of patchy particles. A patch is defined by a solid angle 
with half opening angle S about an axis Ui. Two particles attract each other 
if the angle between two patches in their surfaces is within a given range. 

Our model consists in a continuous version of the previous 6-patch model. 
Figure [8] shows a sketch of the patchy particles. The interparticle potential is 
composed of an isotropic interaction and a directional dependent interaction: 

u{r, n) = eu^^ir) + eX"^''(0/(^) (12) 

1/^^ is the Lennard- Jones potential and the u^~'^^{r) potential is 



i/^-"^(r) 



cr\2a /a\<^ 
r J \r 



(13) 



The angular modulation of the interaction taking into account the alignment 
of the patches is given by 



fm- 



1 9i,9j<6 
otherwise 



where 9i is the angle between any patch and the line joining the center of 
mass of the particles and 6 is the patch half opening angle defining the size 
of the patch. Figure |9] shows the potential we are considering in this work. 
In order to characterize the dynamics of the cluster formation, we will 
determine how S(q,t) changes with time. We will also interested in char- 
acterizing the structure of the clusters that form by calculating the local 
bond-order parameter qQ{i). The quantity is a convenient measure of the 
local crystal order. 
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Figure 9: Plot of the proposed potential. The isotropic part of the potential 
is formed by the Lennard- Jones potential (dashed line). The patch-patch 
interaction is modeled by the a — 2a potential with a — V$> and a cut-off 
of Tc = 1.1 (dotted line). When the patches in two different particles are 
aligned, and the distance between the particles is less than r^ the total 
potential resulting is the combination of the isotropic and the directional 
parts (solid line). 



A preliminary result indicates that the range of the isotropic part of the 
potential plays an important role in the initial states of cluster formation. 
We observed from simulation that the long range isotropic potential [vcut — 
2.5) enhances and induces the particles to aggregate in a short time when 
compared with the short range isotropic potential {vcut — 1-15). Figure 10 
shows the energy per particle and the cluster distribution for both the long 
and short range isotropic parts of the potential. 
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Figure 10: (a) Energy per particle for two different ranges of the isotropic 
part of the potential. In the short range case (open circles), the energy per 
particle slightly fluctuates around the value u/ksT ^ —1.75 while in the 
long range case (filled squares) the energy exhibit larger fluctuations around 
a lower value, (b) Cluster distributions for the two previous conditions. 
The long range potential (right panel) induces the formation of a big cluster 
from the early stages of the simulation. On the contrary, for the short range 
potential (left panel), we only observe a small size cluster distribution. 
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